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The wormlike chain (WLC) model of DNA bending accurately reproduces single-molecule force-extension 
profiles of long (kilobase) chains. These bending statistics over large scales do not, however, establish a 
unique microscopic model for elasticity at the 1-10 bp scale, which holds particular interest in biological 
contexts. Here we examine a class of microscopic models which allow for disruption of base pairing (i.e., a 
'melt' or 'kink', generically an 'excitation') and consequently enhanced local flexibility. We first analyze the 
effect on the excitation free energy of integrating out the spatial degrees of freedom in a wormlike chain. Based 
on this analysis, we present a formulation of these models that ensures consistency with the well-established 
thermodynamics of melting in long chains. Using a new method to calculate cyclization statistics of short 
chains from enhanced-sampling Monte Carlo simulations, we compute J-factors of a meltable wormlike chain 
(MWLC) over a broad range of chain lengths, including very short molecules (30 bp) that have not yet been 
explored experimentally. For chains longer than about 120 bp, including most molecules studied to date in 
the laboratory, we find that melting excitations have little impact on cyclization kinetics. Strong signatures 
of melting, which might be resolved within typical experimental scatter, emerge only for shorter chains. 
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I. INTRODUCTION 

The challenge of constructing coarse-grained models, 
e.g., for large biomolecular systems, lies in accounting for 
the influence of fluctuations that are not explicitly rep- 
resented. Systematic procedures for coarse-graining per- 
form such an accounting, often approximately, so that 
physical consequences of fluctuating solvent densities, 
electric fields, etc. can be incorporated at low computa- 
tional cost. In the simplest cases, integrating out certain 
degrees of freedom just renormalizes interaction parame- 
ters for the remaining variables. In general new types of 
interactions are introduced as well. 

This familiar process of renormalization has impor- 
tant implications for the way experimental data should 
be used to parameterize microscopic models. Measured 
statistics of a particular variable include the effects of all 
other fluctuating degrees of freedom. Thermodynamic 
parameters inferred for that variable do not directly in- 
form the energetics of models that explicitly include other 
degrees of freedom. In some situations these considera- 
tions are transparent: an experimentally derived implicit 
solvent potential should not be used in models that ex- 
plicitly represent solvent fluctuations. In other cases the 
issue can be somewhat more subtle. This paper demon- 
strates the importance of disentangling renormalization 
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effects in the specific context of models for bending of 
nucleic acid molecules. 

The wormlikc chain (WLC) model has proven remark- 
ably successful in reproducing experiments probing the 
conformational flexibility of DNA. In a discretized form 
the model envisions DNA as a fluctuating chain of dis- 
crete links, each inextensible in its length^l, with the 
persistence length £p ~ 50 nm (~150 basepairs [bp]) set- 
ting the contour length scale over which orientational 
correlations decay. The model shows very good agree- 
ment with single-molecule force-extension me asurements 
on kilobase-long A-phage genomic DNApEI The WLC 
with £p = 50 nm also accurately predicts cyclization rates 
of medium-length (hundreds of bp) DNA chains in liga- 
tion experiments!^ In general the WLC model reproduces 
the results of experiments whose observations are domi- 
nated by conformations typical of thermal equilibrium. 

Yet in biological contexts DNA is often bent on much 
shorter length scales than the micron lengths of A-phage 
DNA: prokaryotic transcription initiation^ nucleosomal 
genome compaction in eukaryotespl DNA-binding by ar- 
chitectural proteins such as IHF, 7 and viral DNA pack- 
agingpl all feature DNA bending hundreds of degrees on 
length scales of tens of bp. Furthermore, DNA is in- 
creasingly being used as a prog rammable template for 
constructing nanomaterials, 9 12 where its bending flexi- 
bility over tens of bp can strongly influence the resulting 
structures.^ Knowledge of the mechanical properties of 
short DNA chains is thus essential to understand the role 
of DNA looping in gene regulation, the nature of protein- 
DNA interactions, the pressure generated in viral DNA 
packaging, and the patterned nanomaterial dictated by a 
given DNA scaffold. 

Many microscopic bending potentials produce long 
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length scale bending statistics identical to the WLCP 
and thus the relatively well-established long length scale 
results do not distinguish between substantially different 
possibilities for smaller length scale bending potentials. 
In particular, the energetics associated with large end- 
to-end forces, or equivalently, during exceedingly rare 
large thermal fluctuations, could take on numerous differ- 
ent forms without substantially changing the long length 
scale bending statistics. 

Indeed, recent experiments have suggested that 
oligomers shorter than a single persistence length may 
have considerably different mechanics than a WLC with 
£p = 50 nm. Measuring the end-to-end distance of DNA 
free in solution by Forster Resonance Energy Transfer 
(FRET) and its radius of gyration by small angle X-ray 
scattering (SAXS), Yuan, et al. have inferred an ap- 
parent persistence length of 20 nm for chains comprising 
fewer than 21 bp though it is not clear how this conclu- 
sion is consistent with the accepted value inferred from 
experiments on longer DNA chains. Cyclization rates 
of short chains are particularly dominated by rare large 
fluctuations, and thus should provide a sensitive test of 
DNA's detailed short length scale bending statistics. In- 
deed, DNA ligation experiments by Cloutier and Widom 
found that DNA of approximately 100 bp cyclized as 
much as five orders of magnitude more readily than ex- 
pected from corresponding WLCs with £p — 50 nmP 
By contrast, in similar experiments Du, et al. found the 
cyclization efficiency of small DNA molecules to be in 
good agreement with the traditional WLC model of DNA 
bending albeit with a slightly shorter persistence length 
(£ P , 47 nm)P 

It is easy to imagine that the collective deformations 
of a uniformly buckled rod, when subjected to extreme 
bending, could be superseded by localized excitations 
that render short stretches of the chain very pliable. 
Recent theoretical work has sought to account in this 
manner for the anomalous rates of Cloutier and Widom 
by positing a wormlike chain model incorporating such 
thermally-excited 'melts' or 'kinks'. The apparent per- 
sistence length (after integrating over the excitations) of 
such a model is greater than that of the conventional 
wormlike chain: melting influences bending. In this work 
we demonstrate that, conversely, it follows that bend- 
ing influences melting: the apparent free energy of form- 
ing such an excitation (after integrating over the spatial 
degrees of freedom) is less than the 'bare' free energy 
directly entering into the model Hamiltonian. In §IV| 
we show that as originally parameterized, the meltable 
wormlike chain theory of Yan and MarkcS^Eni improp- 
erly accounts for this entropic stabilization of enhanced- 
flexibility excitations. Specifically, with the original ther- 
modynamic penalty for local melting, that model ap- 
pears to be inconsistent with well-characterized thermo- 
dynamics of bulk DNA melting. We discuss this double- 
counting of certain entropic gains due to melting, and in 
jfV| propose a revised parameterization for the 'meltable' 



uating J-factors that characterize cyclization kinetics, 
from enhanced-sampling Monte Carlo (MC) simulations. 
This method could be applied to any comparably coarse- 
grained model and, unlike existing methods, incurs sim- 
ilar computational cost for both long and short DNA 
chains. In §VIII| we explore cyclization rates using these 
MC simulations, finding that such excitations, properly 
parameterized by well-established DNA melting thermo- 
dynamics, produce no practical difference in cyclization 
rates for DNA chains longer than ^120 bp. Only for 
shorter chains do models including melting excitations 
produce significantly higher cyclization rates than the un- 
mcltable WLC model, though this rate enhancement is 
lower than that seen in some experiments. 



II. THE WORMLIKE CHAIN MODEL WITH THERMAL 
EXCITATIONS OF ENHANCED FLEXIBILITY 

The WLC model envisions DNA as a single chain of 
uniform-length segments (each representing a fixed num- 
ber of basepairs) connected at a set of nodes, where the 
bend angle at a given node i is allowed to fluctuate ac- 
cording to a bending energy quadratic in the local cur- 
vature 
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Here j3 = (fes? 1 ) -1 , £p is the persistence length, tj is the 
unit-length vector pointing from node i — 1 to node i, 
d is the contour length separating adjacent nodes, and 
6i = cos _1 (t i+ i • ii) is the local bending angl e (Fig , [lji). 

The meltable WLC (due to Yan and MarkcPHflf) a dds 
an additional set of variables, the hybridization states of 
those nodes: each node has two internal states, notionally 
corresponding to hybridized and melted local structures. 
Fluctuations in base pairing in this model are represented 
by changes in the hybridization state. The free energy 
associated with a single node in hybridization state m 
and bent at an angle 9 is 
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WLC model. In I VII we outline a new method for eval- 



where ip 1 is the persistence length of hybridized DNA and 
Ip 1 is the persistence length of molten DNA. We have 
defined m = as the hybridized state and m = 1 as the 
molten state. The Kronecker delta, 5i t j, is if i ^ j and 
1 if i = j. The thermodynamic penalty for melting A/i 
is the reversible work required to disrupt base pairing at 
a node of the chain that has a fixed bend angle 9 = 0. 

We also compare the MWLC with two other mod- 
els originally formulated from similar perspectives. The 
kinkable wormlike chain (KWLC)— posits a completely 
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flexible molten state, £p = 0. The spontaneous bend 
model maintains some rigidity in the molten state but 
introduces a preferred bend angle cos -1 7: 
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Nevertheless, the specification of a preferred bend an- 
gle seems less physically plausible for bare DNA than 
for protein-bound DNA, and the assumption of com- 
pletely flexible molten regions is unnecessarily extreme 
for our purposes. Thus we focus primarily on the role 
of thermally-excited local regions of enhanced flexibility 
('melts') in the cyclization kinetics of DNA. 

In these models a chain of contour length (N + l)d has 
N + 2 nodes (which we number to N + 1) connected 
by N + 1 fixed-length links, and N bending degrees of 
freedom (Fig.[l^i). Lacking bending energy, the first and 
last nodes can be ignored in writing the chain's total 



energy: 
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Because nodes contribute to £ c hain in an additive fashion, 
they fluctuate in a statistically independent manner (in 
the absence of collective constraints). 

In considering a model chain that is inextensible along 
its contour, we neglect potential effects of stretching fluc- 
tuations on cyclization. This approximation is supported 
by esti mates of the dsDNA stretch modulus, ^1000 
pN J2122I anc j by the threshold force, ~65 pN, for the over- 
stretching transition to a form ^1.7x longerPBoth val- 
ues well exceed the typical range of forces anticipated in 
the course of loop closure. This model only implicitly 
accounts for structural aspects on length scales smaller 
than a single link, such as the detailed atomic structure 
of the DNA basepair. 



for dimerization fcdim = hm^o [d£>(i)/d£] / Mq, controls 
in certain limits for the ligase enzyme's propensity to 
ligate apposed sticky ends, yielding a measure of chain 
flexibility, the so-called J-factor, J = 2fc cyc /fcdi m - Here 
C(t) and D(t) are the concentrations of ligated circular 
monomers and ligated linear dimers, respectively, and Mq 
is the initial concentration of unligated monomers. 

In the case of low ligase concentration, ligation should 
proceed at a rate proportional to the equilibrium popu- 
lation of precursors adopting ligation-competent confor- 
mations. We define a monomer conformation as ligation- 
competent when one end lies within a critical volume 
v* (of critical radius R*) of the other end of the same 
molecule, while the relative polar and azimuthal angles 
of closure (On — &o and 4>n — 4>o, respectively) lie below 
small threshold values (9* and 0*, respectively). There- 
fore at short times the concentration of ligated circular 
monomers will increase with rate 
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Here k is the ligation rate for ligation-competent 
monomers, and P(|Rat — R | < R*,cos(9n — 9 ) > 
cos 9*, \ 4>n — 0o I < 0*) is the probability that the chain 
is considered ligation-competent. The concentration of 
ligated linear dimers will correspondingly increase as 
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where 



P(\R N - R | < R*,cos(9 N - 60) > cosfT, 

1 — cos 



III. CYCLIZATION KINETICS 

For decades, cyclization reactions have been a popu- 
lar experimental tool for measuring DNA bending elas- 
ticityP A population of DNA chains of given length and 
sequence, with 'sticky' single-stranded overhangs on each 
end, are incubated in solution with a DNA ligase, which 
can covalently link the two ends together when they come 
into spatial proximity and appropriate relative orienta- 
tion. Samples are taken at periodic time intervals and run 
on a gel to quantify the fluorescence from bands corre- 
sponding to populations of circularized monomers, linear 
dimers, and so on. Dividing the rate constant for cycliza- 
tion fc C yc = hm t _j.o [dC(£)/dt] / M by the rate constant 



is the probability of one end of a given molecule lying 
within the critical volume of the end of another molecule, 
within the critical relative polar and azimuthal angles. 
The restriction enzymes used to generate sticky ends in 
the experiments discussed here (Eagl and Hindlll) oper- 
ate on reverse-palindromic sequences, so each monomer 
has an identical overhanging sequence on each end, and 
thus either end can be ligated to another monomer, pro- 
ducing the factor of 2 in Eq. ^ 

We assume that end-to-end distances at which ligation 
occurs are smaller than any length scale characterizing 
chain statistics, and thus that the probability density is 
uniform and isotropic within the critical volume and crit- 
ical angles, 
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(a) I , (b) 




FIG. 1. Geometry of a meltable WLC comprising N + 1 links, (a) Chain cyclization amounts to imposing the constraints 
Ro = Rjv, #o = On, and <f>o = 4>n; where <j>o an d 4>n are azimuthal angles (not shown) corresponding to the polar angles 6q 
and On, respectively. Note that the angle 8o is defined with reference to the vector Rjv — Rjv-i, which is oriented vertically 
in this example configuration, (b) To determine the propensity of apposed chain ends to attain similar relative orientation, 
biased simulations (see j |VI[ | calculate the conditional probability distribution of cos^joia, given Rjv = Ro, where 9j ; n is the 
angle between the vectors Rjv — Rjv_i and Ri — Ro- 

P(|Rat-R | < R*,cos(9 N -9 Q ) > cos0*, \4> N -<j>o\ < <f>*) = v*(1-cob0*)4>* (<5(Rjv-Ro) 5(l-cos[6 N -6 ]) S(^> n -<Po))n 

(8) 



where the angled brackets with subscript N denote a 
canonical average over the 27V angular degrees of free- 
dom (two at each internal node) of an (7V + 2)-node chain 
according to the specified energy function (in our case 
-^chain hi Eq. Q). The Dirac ^-function S(x) is when 
x 7^ and integrates to unity over any region including 
x = 0, thus the latter two (5-functions vanish unless the 
tangent vectors (Ri — Ro)/d and (Rjv+i — RaO/^ are 
identical. (See Fig. [l^,.) 



Under these conditions, the J-factor reduces to an ef- 
fective probability density of the unligated chain forming 
a transient closed loop, with zero relative polar and az- 



r 



imuthal angles: 



J = 2 



2M a lim 
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dC(t) /dD(t) 
dt I dt 
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v*(l — cosd*)(f)* 

x P(\R N - R | < R*,cos(9 N - 6 ) > cos9*, 

|0jv - 0o| < 0") (9c) 

4tt (5{R n - R ) 5(1 - cob[B n - ]) S(0 N - O )) N . 

(9d) 



In {VII we compute the probability density (<5(R 
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Ro) 5(1 — cos[0jv — 6o])5((f)]y — (/)o))n for the meltable 
WLC model outlined in £|TTJ 

When ligase concentration is not sufficiently low, the 
ligation rate of ligation-competent conformations be- 
comes comparable to the rate of their formation from 
ligation-incompetent precursors. In this concentration 
regime ligation rates will thus depend on the dynam- 
ics of formation of ligation-competent conformations, 
not just on their equilibrium probabilities. For circu- 
lar monomers, this reflects the rate of motion in the co- 
ordinate of the end-to-end distance. For linear dimers, 
kinetics will be diffusion-controlled. Thus the ratio of 
circularization and dimerization rates no longer sim- 
ply reflects an equilibrium probability, and the J-factor 
no longer reflects equilibrium free energies. Du, et al. 
showed that at 21°C the threshold ligase concentration 
was ^100 units/ml, and thus argued that Cloutier and 
Widom's ligase concentrations of 150-250 units/ml were 
sufficiently high that they did not measure equilibrium 
bending propensities. Recent experiments by Forties, et 
al.^1 complicate the picture by indicating a threshold at 
37° C of ^400 units/ml, suggesting that Cloutier and 
Widom's experiments at 30° C might have near-threshold 
ligase concentrations. See also Peters and Maher's review 
for further discussionP3 

In the limit of high ligase concentration, the rate- 
limiting step is formation of a ligation-competent confor- 
mation, and hence the J-factor would reflect the effective 
reaction rate of aligning the ends. This quantity can in 
principle be calculated from molecular dynamics simu- 
lations of models with more detailed representation of 
the DNA: though the relevant dynamical modes might 
be complex, a first approximation might assume diffu- 
sive motion in the end-to-end distance and join angle, 
and subsequently solve for the diffusion-controlled steady 
state using a Langevin or Fokker-Planck approach. Com- 
parison with experiments at high ligase concentration 
could thus yield insights into DNA dynamics on these 
short length scales. Cheng, et al. recently conducted 
similar experimental and computational examination of 
the cyclization dynamics of single-stranded DNAP-3 In 
the rest of this paper, we restrict our discussion to the 
regime of low ligase concentration. 



IV. CONNECTION WITH LATTICE MODELS 



The chain's total free energy in this description sim- 
ilarly depends only on the hybridization variables, 
-EchainljwJ] = J2^Ei m i)- Implicit in the free energy 
assignment of Eq. ( |10[ ) is an integration over conforma- 
tional fluctuations weighted by the Boltzmann distribu- 
tion, 



exp(-/3E cha i n [{ni}]) 



(11) 



d(cos^i) / <i(cos#2) ■•■ / d(cos0/v-i) 
l J-i J-i 

x exp(-PE cha . in [{ni},{0i}]) . 



Analytical evaluation of this integral is straightforward 
for the meltable WLC, since the Boltzmann weight fac- 
torizes. 

The partition function at a given node for hybridiza- 
tion state to (integrating Eq. ^ over all bend angles 
cos 9) is 
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^ + 5 m , 1 (^-l|) eXp{ ~ ,5 - l/3AM} (12b) 
'l-expj-^ + ^x^-^)] 

The free energy difference between melted (to = 1) and 
non-melted (to = 0) states is AF = F(l) ~ F(0) = 
—k^T ln[g(l)/g(0)]. We arrive at a relationship between 
the 'bare' free energy of melting (A/i MWLC ) and its 
'renormalized' counterpart (A/i) that accounts for the 
influence of bending fluctuations: 
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The second term in the RHS of Eq. (13a I accounts 



for the greater range of bending motion available to a 
molten segment of the chain, which is considerable. For 
the length parameters employed by Yan and Marked 
(£p = 50 nm, if = 1 nm, d = 1 nm), the renormal- 
ization evaluates to 



One-dimensional lattice models for fluctuations in 
DNA hybridization have a long history in the study 
of bulk melting thermodynamics.^ Systems of interest 
are typically free of any external mechanical constraints, 
so the models typically make no statement about chain 
structure. In particular they assign free energy to a short 
segment of the chain based only on its hybridization state 
nr. 
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Am- A/x MWLC « -3.8fc B T 
Similar analysis for the KWLC yields 
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and for the spontaneous bend model 
0Aji = (3A/j, SB - In 
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Conversely, integrating over mel t deg rees of freedom 
renormalizes the persistence lengthPSGSl 



To maintain a given apparent persistence length £p 
for a given renormalized melt free energy /3A[i, the hy- 
bridized persistence length £p is thus a function of the 
melt persistence length 



«H,MWLC 



For the KWLC model we find 
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V. PARAMETERIZATION 

The large body of work on DNA melting provides good 
estimates for the thermodynamic parameters appropriate 
to one-dimensional lattice models such as Eq. (10). Bulk 



melting free energies include a renormalization [Eq. ( 13 )] 
due to integrating out bending fluctuations, therefore 
proper use of melting data to parameterize microscopic 
models must involve disentangling this renormalization. 
If each node or lattice site represents three basepairs, 
then we expect A/x « 8 — llk B T at ambient condi- 
tions.^ Yan and Marko presented an argument for this 
range of values,^ but no distinction was made there 
between A/i MWLC and A/i. If we assign the inferred 
melting penalty to the bare free energy rather than the 
renormalized value (i.e., A/i MWLC w 8 — 11 k^T and 
A/i « 4 — 7 k-gT), hybridization statistics of the meltable 
WLC at room temperature are not consistent with ex- 
periment. Specifically, a three-base region is melted with 
a probability between 7.2 x 10~ 4 and 1.4 x 10~ 2 (i.e., 
between e -4 2 and e~ 7 - 2 ), instead of between 1.7 x 10~ 5 
and 3.3 x 10~ 4 (i.e., between e -8 and e -11 ) as suggested 
by thermodynamic measurements. 

Furthermore, this significantly-increased probability of 
melts yields chains that are too compliant to be consis- 
tent with longer length scale force-extension experiments. 
The persistence length is determined by (cos 8), the av- 
erage cosine of the angle at a given node. For a meltable 
node this average cosine is an equilibrium-weighted lin- 
ear combination of the average cosines of a melted or 
hybridized node, 



1 



(20) 

producing for the MWLC an apparent persistence length 
(averaged over melt fluctuations) 
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In particular, for the range of bare melt free energies 
A/z MWLC = 8 — 11/cbT, the persistence length is renor- 
malized to ^mwlc _ 3 q_ 48 nirL 



and for the spontaneous bend model with at least modest 
melt persistence length 
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Based on Eqs. ( 14 



alistic value for the 



and ( 22 ) we determine a more re- 



Dare free energy of local melting in 
the MWLC, A/i MWLC « 12 - 15 k B T, which is consis- 
tent with bulk thermodynamics by construction. The 
increased A/z MWLC , designed to avoid double-counting 
entropic gains due to increased flexibility, will clearly act 
to suppress melting and to diminish its importance in the 
kinetics of cyclization. We have quantified this suppres- 
sion using computer simulations. 



VI. RELATING /-FACTORS TO DISTRIBUTIONS OF 
END-TO-END DISTANCE AND ANGLE IN A 
DISCRETIZED CHAIN 



The J-factor was expressed in §III| in terms of an equi- 
librium probability density that lends itself to straight- 
forward evaluation in computer simulations. For a chain 
of contour length (N + l)d, whose N + 2 nodes are lo- 
catedat positions {Ro, Hi, • • • , R-jv+i}) we rewrite 
Eq. (l9d| as 



J(N) = 4tt- 



i(loop) 
!_N 
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(25) 



We have defined partition functions Q^ r °° p ' > for a cyclized 

(N + 2)-node meltable WLC, and Q$ = q N for the 
same molecule absent end constraints. The single-node 
partition function q can be evaluated as 



q = 2n d(cos9) e 



-/3E(cos0) 



(26) 



for apparent bending energy 



E(cos6) = -k B T \ n J2e~ 0E[m ' cose) - (27) 

n=0 
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Thus, 



VII. MONTE CARLO SIMULATIONS 



J(N) = 4tt / d(cos 6»i) " " " y d ( cos °n) (28) 
x / d</>i ••• / d</>Ar- 



Eq. ( 25 ) , together with an expansion of the Boltzmann 
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x ^(Rjv - R ) 6(1 - cos[6 N - 9 \) 5(<j> N - fa) 

The first 5-function imposes no constraint on the vari- 
ables cos On or <f>N (see Fig. [jj,), and thus integrating 
over both is trivial, giving 



J(N) = iwq' 1 ( S(K N - R ) e -^(c°.9o)\ . (29) 
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We then use the delta-function identity 

e -0E(cose o ) = (30) 
d(cos 9 join ) e - p&{cos e ^ S{cos e ^~ cos t 

where the join angle #j in describes the relative orienta- 
tions of the chain's ends, cos#j i n = (Rv — RiV-i)(Ri — 
Ro)/cZ 2 , when it is poised for cyclization, Rjy = Ro 
(i.e., R = 0). (See Fig. [T]d.) Finally, we multiply and 
divide by (5(Rjv — Ro))jv-i; substitute the probability 
density at zero end-to-end extension, pjy(Ro = Rat) = 
(S(Rn — Ro)) jv-i, an d introduce the conditional proba- 
bility density at a given bend angle of a particular node 
(given a cyclized chain), 

Mcos0 join |R = O) = (31) 
(S(R N -R )<S(cos# join - cos0 o ))iv_i 
/ (5(R N - R ))jv-i • 

These manipulations enable us to express the J-factor in 
terms of quantities readily estimated from sampling of 
canonically distributed microstates, 

J(N) = Anq- 1 PN (R = 0)J d(cos # join ) (32) 
x e-^( cose -'») PA ,(cos0 join | R = 0) . 



distribution in spherical harmonics, was used by Yan 
and MarkcEO to evaluate J(N) for the meltable WLC. 
Their approach becomes quite numerically demanding, 
however, for the short chains of interest here. We have 
instead used Monte Carlo simulations to estimate the 
probability distributions appearing in Eq. ( 32 ) . This ap- 



proach also permits the simple extension of the model 
to more complicated energetics, involving considerations 
such as volume exclusion and DNA helicityP^ 

We employ two kinds of trial moves to efficiently nav- 
igate conformational space while maintaining the geo- 
metric constraints of the molecule. 'Free' rotations di- 
vide the chain in two by selecting a randomly chosen 
node, then rotate the two segments relative to one an- 
other about a randomly chosen axis through the dividing 
node. Crankshaft moves rotate a segment of the chain 
between two randomly chosen nodes, about the axis con- 
necting them. The number of contiguous nodes rotated 
in crankshaft moves is chosen from a uniform distribution 
on the interval [1,10]. Angles for all rotations are selected 
at random from Gaussian distributions whose variances 
were adjusted until roughly 1/2 of all trial moves are 
accepted. Melt moves attempt to change the hybridiza- 
tion state of a randomly chosen basepair, either from 
hybridized to melted or vice versa. 

Calculating cyclization probabilities for short DNA 
chains via standard equilibrium sampling is computation- 
ally challenging because of the extreme energetic cost of 
juxtaposing the two ends with the correct orientation. 
We estimated the equilibrium average on the right-hand 
side of Eq. ( 25 ) by sampling from the Boltzmann dis- 
tribution using Wang-Landau sampling)22l an adaptive 
method for broadly sampling fluctuations of a chosen or- 
der parameter without any prior knowledge about the 
details or even the gross shape of the probability distri- 
bution. 

This biased sampling approach allowed us to construct 
the distribution p^(R) = (S(H — (Rjv — Ro)))iV of the 
end-to-end vector R for an uncyclized chain, over the 
full range of |R| (Fig. [2]). Similarly biased sampling, this 
time of a cyclized chain with no bending energy at one 
node, was used to construct the conditional angular dis- 
tribution at that freely-jointed node (Fig. [3]): 



Pat (cos ( 



R = 0) 



(S (cosfljoin - cr 2 [R! - R ] • [Rjv - Rjv-i]) S(R N - R )) 



(5(R N - R )) 



JY 



(33) 



For each chain length and melting free energy, 3 different 
sampling trajectories were run from identical initial con- 
ditions but with different random number sequences to 
calculate uncertainties. Each Wang-Landau iteration ran 



until each of the 1000 spatial bins (or 200 angular bins) 
was visited at least 100 times. The 10 spatial bins at 
greatest extension were excluded from this criterion due 
to their exceedingly small probabilities. The increment 



to the Wang-Landau bias began at 1 k B T and decreased 
by 10% at each subsequent iteration, until the increment 
reached 10~ 7 k B T (10~ 8 k B T for angular simulations of 
chains shorter than 120 bp). A few angular simulations 
of shorter chains achieved greater precision with a 20% 
(rather than 10%) decrease in increment per iteration. 

Probability distributions of the end-to-end distance R 
and join angle 0j o i n , for meltable wormlikc chains of vary- 
ing length and melt free energy, are shown in Figs. [2] 
and [3] For meltable chains, the end-to-end distribution 
approximates at large distances that of the un-meltable 
chain. As end-to-end distance decreases below the range 
of typical fluctuations of a WLC, melting begins to sig- 
nificantly facilitate achieving loop closure of sufficiently 
short chains. Their significant energetic cost requires 
that substantial bending forces develop before such lo- 
cal excitations become favorable. 

The general shape of pjv(R), which has been discussed 
at length elsewhere,^ can be simply understood in terms 
of entropic elasticity at large R and bending energy at 
small R. Straightening a chain imposes high entropic 
cost, so that probability decreases sharply with R at large 
extension. Imposing a short end-to-end distance requires 
significant bending and thus high enthalpic cost, explain- 
ing the low probability for very bent chains. The longest 
chain is easiest to cyclize since the total bend can be 
distributed over the largest number of basepairs. 

At the very shortest end-to-end distances, meltable 
chains show a slight increase in pjy(R) as distance de- 
creases, for geometric reasons. Imposing such a tight 
bend essentially requires a melt at the middle node of the 
chain; once that node is melted it serves as a nearly free 
joint between two nearly rigid chain segments on either 
side. Taking the joint as the origin of a molecular ref- 
erence frame, the two endpoints thus effectively explore 
the surface of a sphere. As they approach one another in 
this two-dimensional geometry, the number of accessible 
conformations with a given value of R grows only lin- 
early, n ~ R (rather than quadratically, as it would for 
a very flexible chain). The probability per unit volume 
Pat(R) thus scales here roughly as fi(i?) / '(AirR 2 ) ~ l/R. 
Because the chain ends are not strictly confined to a two- 
dimensional surface, pjv(R) does not diverge as R — » 0, 
but it does grow over a small range of R. A similar sce- 
nario, now involving melts at nodes one or two removed 
from the central one, accounts for the small oscillations 
at slightly longer end separations. Consistent with this 
explanation, the local maxima in probability appear at 
end-to-end distances that are even integer multiples of 
the link length. 

For unmeltable chains, the join angle statistics (for a 
chain that meets at its ends but lacks bending energy at 
that meeting point) in Fig. [3] show a probability max- 
imum at intermediate angles. This stems from the en- 
thalpic preference for equally distributing the bend over 
the entire contour of the chain, producing a teardrop- 
shaped chain contour with join angle intermediate be- 
tween 0° (straight) and 180° (fully bent) P Meltable short 



chains differ substantially. In accommodating loop clo- 
sure, the energetic cost of severe bending distributed 
across the entire uniformly curved chain (without any 
melts) is substantially higher than that of a single melt 
at the middle basepair opposite the junction. The lat- 
ter state is energetically minimized with a join angle of 
^180° and very little bending on either side. Conversely, 
the low probability of a join angle of O°(cos0j O in = 1) re- 
flects the need for formation of a second melt, essentially 
required to bring the ends of such short chains into com- 
mon orientation. 



VIII. J-FACTORS 

We have confirmed numerically that this MC approach 
is consistent with the spherical harmonic expansion used 
by Yan and Marko. 1 ^ For the chains we compared, of 
length 135-225 bp, the deviations between their results 
and our calculations are not statistically significant (data 
not shown). Henceforth we calculate MWLC J- factors 
for parameters that differ in two respects from those of 
Yan and Marked (i) melted sections are slightly stiffer, 
^p 1 = 2.5 nm, than in their calculations (where ip = 1 
nm); and (ii) the range of melting free energy values we 
explore extends to the stiffer thermodynamic penalties 
suggested in SjV] We argue that the longer persistence 
length of unhybridized chains is more consistent with ex- 
perimental estim ates for single-stranded DNA (ssDNA) 
bending elasticityP^U since a melted section composed of 
two single strands of DNA in close apposition would be 
expected to have a persistence length at least twice that 
of ssDNA. 

All of the model results we report, like those of Yan 
and Marko, neglect the helical nature of hybridized ds- 
DNA, implicitly assuming that two ends juxtaposed with 
any helical phasing are ligation-compctcnt. The calcula- 
tions thus represent an upper bound for the J-factor of a 
helical wormlike chain polymer, a bound that is tightest 
for chains that are an integral number of helical repeats 
in length. 

Fig. [4] shows our numerical results for contour lengths 
of 10-120 nm, with melting energies from Yan and 
Marko 19 (a) or with values of A/i that are consistent 
with bulk melting thermodynamics (b). The shift from 
A/i = 8 - llk B T (Fig. fy) to Afj, = 8 - llfe B T 
(Fig. |4Jd) represents an increase in melting free energy 
of 3.7 k B T and hence considerably suppresses melts. In 
fact, whereas for A/i = 8 — 11 k B T only the 360 bp chain 
J-factors are unaffected by melts, for A/i = 8 — 11 k B T 
(A/x = 11.7 — 14.7 k B T) meltable chains as short as 150 
bp cyclize with propensities indistinguishable from the 
melt-frcc limit, A/i = A/i = oo. For shorter chains the 
meltable models produce a J-factor significantly greater 
than that for the melt-free (A/i = A/i = oo) model, 
because cyclization dictates a bend so sharp that the en- 
ergetic gain from localizing the bend at a floppy melted 
section overcomes the energetic cost of forming the melt. 



9 




20 40 60 80 100 120 10 20 30 



R, nm R, nm 

FIG. 2. Bending statistics for meltable wormlike chains of varying contour length and melting thermodynamics. Error bars 
show ± 2 standard errors of the mean (s.e.m.), calculated from three independent simulations, (a) Distributions of end-to-end 
distance R, for chains comprising 30 — 360 bp, generated by Wang-Landau sampling. Melts are disallowed here, i.e., /3A/i = oo. 
(b) Distributions of end-to-end distance R, for 90 bp meltable wormlike chains for several values of the (renormalized) free 
energy of melting, Afi. 



Surprisingly, at the short chain lengths (less than 
~90 bp), a meltable chain becomes easier to cyclize 
with shorter chain length, reflected in a higher J-factor. 
This trend is rationalized through consideration of the 
enthalpic and entropic costs of cyclization. Cyclizing 
a short fully-hybridized DNA is so energetically costly 
(seen in the sharp drop in J-factor for the un-meltable 
chain (black curve)) that a short meltable chain essen- 
tially requires a melt near the middle to cyclize. Such 
a cyclized structure with a single melt will have an en- 
thalpy largely independent of chain length, due to the 
high compliance of the melted region. However the en- 
tropic cost of cyclization decreases with decreasing chain 
length, so for short meltable chains the probability of 
cyclization increases with decreasing chain length. 

This striking increase in J-factor for very short chains 
also occurs for the KWLCP It does not appear in the 
spontaneous bend model with sufficiently rigid molten 
sections, nor in a model that fixes a particular bend an- 



gle^ (equivalent to the spontaneous bend model with 
infinite molten persistence length £p), because in these 
models cyclization of shorter chains still requires tighter 
bending and thus greater enthalpic cost. It also would 
not occur for the linear subelastic chairJ^ as that model 
does not concentrate bending at one location. No ex- 
periments have yet probed the short lengths where this 
phenomenon is seen for the MWLC or KWLC. 

We now compare with experiment these computational 
results for meltable wormlike chains with the stiffer renor- 
malized melting free energies A/i = 8 — 11 fee J (Fig. [4Jd 
inset) . Cloutier and Widom's experimental cyclization 
ratea^ show significant scatter, presumably reflecting 
sequence heterogeneity in melt- and bend- propensities 
even among 'random' sequences (which exclude nucleo- 
some positioning sequences known to have unusual bend- 
ing characteristics). The meltable wormlike chain model 
predicts cyclization rates significantly lower than those 
inferred from experiment over the range of 90-120 bp, 
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FIG. 3. Statistics of chain's initial orientation (60), relative to that of its final link (On), for cyclized meltable wormlike 
chains of varying contour length and melting thermodynamics. Error bars show ± 2 s.e.m. (a) Distributions of cos0j o in, for 
wormlike chains comprising 30 — 360 bp, generated by Wang-Landau sampling. Melts are disallowed here, i.e., /3Afi = 00. (b) 
Distributions of cos(9join for 90 bp meltable wormlike chains for several values of the (renormalized) free energy of melting, A[i. 



which includes the majority of Cloutier and Widom's 
sequences. Melting excitations, when made thermody- 
namically consistent with bulk melting behavior, appear 
insufficient to explain their results. 

The experimental results of Du, et al. 17 highlight 
the J-factor's strong helical dependence at these short 
lengths. Local maxima in J occur for chain lengths that 
are integer multiples of the helical repeat, thus requiring 
minimal over- or under-twisting to bring apposing ends 
into helical register. These chain lengths should offer 
the most straightforward comparison with the meltable 
wormlike chain model (which implicitly assumes perfect 
helical register). Throughout the range of N studied in 
these experiments, these maxima are slightly higher than 
corresponding J-factors of the meltable wormlikc chain, 
or compare favorably with results for the smallest val- 
ues of A/i. Given debate over the precise value for DNA 
persistence length within the range of ^45-53 nm, the 
experiments of Du, et al. appear to be consistent with 
a meltable wormlike chain model. But in this range of 



uncertainty, their results are almost equally consistent 
with an unmeltable wormlike chain, which has similar 
J-factors to the meltable wormlike chain for molecules 
longer than ~120 bp. Shorter chains would offer a much 
more discriminating test: for N ^ 50 bp, the predicted 
J-factors of meltable and un-meltable chains differ enor- 
mously. Experiments on such molecules should unam- 
biguously address the importance of thermally-excited 
melts in the severe bending of dsDNA. 



IX. TEMPERATURE DEPENDENCE 

Enhanced cyclization of long DNA strands (thou- 
sands of bp) due to elevated temperature was estab- 
lished decades agoF^ Yet there are few experiments (to 
our knowledge) probing temperature- variation of cycliza- 
tion kinetics at the short lengths of interest here. Given 
the expected strong Arrhenius temperature-dependence 
of melting, such measurements should clarify the role 
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FIG. 4. Cyclization kinetics calculated for meltable wormlike chains of varying contour lengths and melting thermodynamics. 
Error bars show ± 2 s.e.m. (a) J-factor plotted as a function of contour length, for bare melting free energies A/i = 8 ksT 
(magenta), 9ksT (blue), WIcbT (green), llksT (red), and in the absence of melts, Afj, — oo (black), (b) J-factor plotted as a 
function of contour length, with renormahzed melting free energies A[i = 8 fcaT (magenta), 9 ksT (blue), 10 ksT (green), 11 ksT 
(red), and in the absence of melts, A/i = oo (black). Solid black line shows predictions from the approximate analytic theory of 
Shimada and Yamakawa (in the absence of melts) Black squares show experimental results of Cloutier and Widomp^gray 
squares those of Du, et a/P-^ 



of dehybridization or even explain discrepancies between 
different experiments. Here we scrutinize these possibil- 
ities through the thermal sensitivity of an MWLC. 

In this section we present numerical results for J- 
factors as a function of T. We assume (i) temperature- 
independent bending rigidities ksT £p and fc B T^p ; and 
(ii) a constant value of the apparent bending rigidity 
kftT £p, a weighted combination of hybridized and melted 
bending rigidities that is dominated by the former. Tran- 
sient electric birefringence experiments in the range of 
20-43° support this latter assumption. 

Fig. [H] shows computed J-factors, for the case of Afi = 
lOfceT, at several different temperatures. We find that 
a change of 10° C has very little impact on J for all 
chain lengths examined. In fact varying temperature 
by as much as 30° C never altered J by more than a 
factor of three, for all lengths examined. Experiments 
spanning 23-42° C found similar temperature-sensitivity 



for a 200 bp segment of A DNA and a 116 bp chain 
with sequence chosen to minimize melting, but some- 
what greater temperature-sensitivity in a 116 bp chain 
with sequence designed to accentuate meltingPS Since 
the temperatures in the respective experiments of Du et 
al. and Cloutier and Widom differ by only nine degrees 
C, temperature variation is unlikely to explain differences 
in their findings. 



X. CONCLUSION 

In this work we have explored the mutual influence 
of distinct degrees of freedom characterizing the micro- 
scopic state of a chain molecule. Our results highlight 
how the presence or absence of one fluctuating variable 
shapes the effective forces on another. From a practi- 
cal perspective, these results caution that grafting an 
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FIG. 5. Temperature dependence of cyclization kinetics for 
meltable wormlike chains. J-factors are plotted as functions 
of contour length for T = 15, 25, 35, and 45°C. A moderate 
melt enthalpy, A/i = 10 ksT, was used in all cases. Error bars 
show ± 2 s.e.m. Solid black line shows predictions from the 
approximate analytic theory of Shimada and Yamakawa (in 
the absence of melts) 

additional degree of freedom onto a model changes the 
apparent parameters governing the fluctuations of other 
degrees of freedom. 

These concepts find a concrete example in models of 
DNA flexibility that incorporate regions of thermally- 
enhanced flexibility It seems immediately clear that in- 
troducing the possibility of flexible excitations to a poly- 
mer will reduce its apparent persistence length, hence the 
invocation of melts in attempts to understand anoma- 
lously high cyclization rates. Somewhat more subtly, 
adding bending degrees of freedom (that differ in their 
natural fluctuations between hybridized and melted sec- 
tions of the chain) to each site in a lattice model for 
melting can substantially renormalize the free energy of 
dehybridization. For the particular parameters used in 
our simulations (£p = 50 nm, £^ = 2.5 nm, d = 1 nm), 
the bending-induced increase in apparent melt free en- 
ergy amounts to 3.7 k^T. 

Based on this distinction between bare and appar- 
ent melt free energies, we develop a thermodynami- 
cally consistent treatment of the meltable wormlike chain 
model that properly accounts for the greater entropy of 
more flexible melted sections, and thus permits accurate 
parameterization from melting experiments. For suffi- 
ciently short DNA chains, the conformations that lead 
to cyclization are exceedingly rare and thus their equilib- 
rium probability is greatly facilitated by the transiently 
enhanced flexibility of thermally-excited melts, even at 
the energetic cost of disrupting DNA basepairing and 
stacking interactions. However, the proper parameter- 
ization of melting energies in the range A/j, = 8—11 k^T 
pushes the onset of such melt-enhanced cyclization to 
significantly shorter lengths, so much so that the cy- 



clization propensity of a meltable chain is practically 
indistinguishable from an un-meltable chain for chain 
lengths greater than ^120 bp. Predictions for meltable 
and un-meltable chains diverge significantly around ^80 
bp, suggesting that experiments on even shorter chains 
would distinguish between models. Indeed at the short- 
est lengths, meltable chains even show an increase in 
cyclization propensity with decreasing chain length, due 
to entropic effects. 

We have also demonstrated that coarse-grained model- 
ing with Wang-Landau sampling can be a useful method- 
ology for calculating J-factors for short DNA chains. 
This basic approach could find further application in ex- 
ploring additional aspects of DNA bending elasticity, of 
great relevance for DNA compaction and protein-DNA 
interactions. Elaborations of the MWLC model may 
be required for these purposes. Our existing compu- 
tational framework could easily incorporate the depen- 
dence of bending elasticity on salt concentrations and 
DNA sequence, through altered free energies of melting^ 
and al tered flexibility of a given melted or hybridized 
sectionP3in] Other possible straightforward amendments 
include more detailed melt thermodynamics (for exam- 
ple incorporating stacking interactions,^ a specific heat 
terrrP^l or spatial and dynamical informational) ; volume 
exclusion, and explicit representation of DNA helicityP^ 
Nevertheless, the meltable WLC still provides a parsi- 
monious model lending important insights regarding the 
effect on DNA cyclization behavior of greater flexibility 
conferred by thermally-excited defects. 
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